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Abstract 

The lattice Bhatnagar-Gross-Krook modeled fluid has an 
unchangeable unit Prandtl number. A simple method is 
introduced in this letter to formulate a flexible Prandtl 
number for the modeled fluid. The effectiveness was 
demonstrated by numerical simulations of the Couette 
flow. 



Lattice Boltzmann (LB) method^, which is a semi- 
ii. continuous version of lattice gas automata (LGA)^, has 
Jheen accepted recently as an alternative numerical tool 
for the simulation of transport phenomena in the fluid 
^system. The latest modiflcation^^® of the LB method 



idoes remove many disadvantages of the anecstral LGA 
^model, such as the lack of Galilean in variance, the veloc- 
^ity dependent pressure, the high noise-signal ratio, the 
high viscosity and the inflexible sound speed, etc. On 
l^the other hand, the inheritated particle nature keeps the 
^^high efficiency for the LB method when applied in the 
^^simulation of flows under complicated geometries and re- 
T-Hserves it as one of the most appealing candidates for the 
^Dcomputing models on the massively parallel frame in the 
'^next decade. 

Selecting a regular spatial lattice and rendering each 
C/^link vector of this lattice into the vector for the parti- 
^l^le flight velocity, a phase space for the particle motion 
I can be deflned. The single particle distribution may be 
04enoted, in the discrete spatial and velocity spaces, as 
QNpi(x,t), where d t are spatial and temporal co- 

Qordinates. Indices, p and i, are for sub-lattices and lat- 
Otice links respectively. The lattice Boltzmann equation 
(LBE) to be solved in the process of the LB simulation 
can be written as. 



Npi(x + Cpi,t + 1) - Npi(x,t) = Ap 



(1) 



Here, Cpi is the link vector, which represents the flight ve- 
locity of particle in the meantime. The Ihs of (1) stands 
for the propagation of particle, while Api is deflned as 
the collision term, which is generically linearized as a 
product of a collision matrix and the perturbative part 
of particle distribution, i.e. {Npi — N^f ). The structure 
and elements of the collision matrix shall certainly en- 
sure the conservation of densities of mass (p), momentum 
(pu) and thermal energy (e). More often, the famous 



Bhatnagar-Gross-Krook (BGK) or the single time relax- 
ation approximation (STRA) is introduced in (1), which 
transforms the collision term as simple as the following 
expression. 
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(2) 



Here r is the relaxational time period, through which the 
particle distribution approaches its equilibrium value. 
When the equilibrium particle distribution is expanded 
in the low speed limit as 



Npi = Ap(p,e) + Mp(p,e)(uaCpia) 

+ Gp(p, e)M^ -I- Jp(p, e)(UaCpiaf + O(u^) . 



(3) 



and a set of constraints for the expansion parameters 
Ap,Mp...Jp is properly chosen, a perturbative multi- 
scale analysis would lead to macro-dynamic equations 
in the compressible Navier-Stokes form, provided that 
the lattice is symmetric enough to ensure the isotropy 
for the flux tensors appearing in these equations. This 
serves as a theoretical proof for the use of the LB or 
the lattice BGK fluid dynamic model and a justiflcation 
for naming such a semi-continuous, interacting system 
the lattice BGK fluid, as in the title of this paper. The 
state equation of the current fluid is identical to that 
of perfect gas, because there is no inter-particle force 
considered. Hence the temperature is proportional to the 
thermal energy, speciflcally in two dimensions, one has 
T = e. Transport properties can be described once the 
shear viscosity and the heat conductivity are identifled. 
Indeed, in the D dimensional space they are written as 
follows , 
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(4) 



(5) 



As in usual hydrodynamic text, i> is used for the shear 
viscosity and k the heat conductivity. The ratio of these 
two quantities will be independent of any thermody- 
namic variable except the space dimension, from both 
(4) and (5). Therefore, the Prandtl number of the lattice 
BGK fluid deflned as ^— , where Cp denotes the speciflc 
heat, would have a sole and unit numerical value. 

The reason for the unit Prandtl number of the lattice 
BGK fluid may be understood through a more intuitive 
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way: the form of the colhsion term in this model (2) in- 
dicates an incorporation of the diagonal matrix —-Sij, 
where 6ij is the Kronecker tensor. Eigenvalues for the 
shear and energy transport modes are equal in this case 
and therefore result in the unit value for the Prantl num- 
ber. 

The frozen unit-value Prandtl number for the lattice 
BGK fluid is not favorite, in sense of the flexibility of a 
fluid model, as in many cases one is interested in observ- 
ing the change in the transport of momentum and heat, 
along with the variation of the Prandtl number. For 
example, a bifurcation from the steady roll solution to 
the time dependent thermal plumes is found in the high 
Rayleigh number flow in the natural convection, when 
the Prandtl number of the fluid switches from a lower to 
a higher numerical value. 

To solve this difficulty in the lattice BGK model, an 
additional parameter a can be included in the aforemen- 
tioned constraint for the velocity expansion parameter, 
specifically for Ap in the following form. 
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(6) 



With a standard derivation, (6) will make the ratio of 
heat and momentum transfer tunable, which can result 
in a (T parameterized Prandtl number as follows. 
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(7) 



To verify this formula, numerical measurements of the 
Prandtl number have to be carried out. In this paper, 
the measurements were realized by simulating Couette 
fiows between two parallel fiat walls of which one is set 
at rest, the other moving with a constant velocity U in 
its own plane. It is also postulated that the temperature 
is constant along the wall, but different for the two walls. 



y = 0:T = To; y = h : T = Ti 



(8) 



Here h is the distance between the two walls and in this 
case temperatures are set to be To < Ti . When the 
fiow reaches the steady stage, the temperature profile is 
measured and compared with the theorectical solution: 



Ti-To 
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(9) 



r]+ -(Pr X E)r?(l - r?) , 



where r? = | is the normalized width for the channel and 
E is the dimensionaless Eckert number defined as 
For Ti = Tq case, the solution is obtained as 
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(10) 



which has a parabolic profile across the channel. 

The results are shown in Fig. 1. The quality of the 
modification in (6) is estimated by inspecting the deva- 
tion between the meaured results and analytically pre- 
dicted cureve. It is obvious from the figure that the cur- 
rent method is only reliable at the neighbourhood a = I. 
In other words, the Prandtl number can only be varied 
little from the coherent unit value. The reason for this 
failure can be disclosed by a further study of the heat 
fiux tensor, which, because of the incorporation of a, 
changed into the following form. 



4(T> + 2) 
T>2 



(11) 



[{2a - l)pedae + (1 - a)e^dap] " ^ • 

The existence of daP is unphysical and causes the dis- 
agreement stated above. In order to do away with the 
effect of the density gradient term, another constraint 
involving a may be considered as follows. 



4:(D+2) 
T>2 



p[{l - a)e^ + ae] 



(12) 



where e is the averaged termperature. The unphysical 
part of the heat fiux tensor then becomes 



(1 - '^)(e' - e') 



daP, 



(13) 



so that it has the minimum effect when the distribution 
of thermal energy in the fiow is uniform enough, namely 
in the small temperature difference limit, provided that 
a is not too far from the unit value. 

The effectiveness of such a formulation can be demon- 
strated by the same numerical measurement of the 
Prandtl number described before. Again results are plot- 
ted in Fig. 1 with different symbols for different temper- 
ature differences. The agreement with (7) is extremely 
well, when a is selected from 0.65 ~ 1.35. Correspond- 
ingly the Prandtl number may vary from 3.3 to 0.59. 
The breakdown of this method had been also observed 
either when a exceeds the range, or when the tempera- 
ture difference is too large, e.g. Ae > 0.1. These may 
be considered as restrictions when the new lattice BGK 
model with variable Prandtl number is to be applied in 
any thermohydraulic problems. 

First application of this model is chosen as the Couette 
fiow used above as the numerical experiment to measure 
the Prandtl number of the modeled fiuid. Temperature 
profiles in steady state are shown in Fig. 2. The Eckert 
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numbers for these simulations are the same, i.e. E = 2.0, 
and the temperature difference on the two walls is set 
to 0.01, keeping the method under the valid condition. 
Subsequently different values for a are provided to vary 
the Prandtl number of the fluid as 0.5, 1.0, 2.0. As the 
product Pr X E changes, the heat flux from the upper 
wall to the fluid changes sign. From the temperature 
profiles in Fig. 2, it is clear that cooling of the upper 
wall occurs in case of Pr x E < 2, heating of the upper 
wall in case of Pr x E > 2, and no heat fiows between 
the upper wall and the lattice BGK fiuid as Pr x E = 
2. Furthermore, numerical results and the Naver-Stokes 
solution are found to remarkably agree with each other 
in the same figure. 

Change of the Prandtl number of a fiuid system in 
quantitative sense by microscopic approaches is nontriv- 
ial, as even for the full and continuous Boltzmann equa- 
tion, it is difficult to find the corresponding coUisional 
cross sections®. In the case of discrete Boltzmann equa- 
tion, collision matrices which have different eigenvalues 
controlling transport modes for heat and momentum can 
be used^ to solve the difficulty. Alternatively, terms in- 
cluding gradients of velocity and temperature may be 
incorporated into the equilibrium particle distribution, 
which is analogous to the formulation of Liu's modified 
countinous BGK equation®. Both methods can be ap- 
plied to the lattice BGK model, though both of which 
have difficulties. First the use of collision matrix with 
different eigenvalues may lose the physically correct form 
for the viscous work in the macroscopic equation for 
the energy conservation^'^. On the other hand, gradi- 
ent terms in the equilibrium particle distribution usually 
lower the efficiency of lattice BGK computation, because 
the absolute locality is destroyed in that information on 
the next nearest neighbours are required. Howevere, the 
method introduced in this paper is both simple and ef- 
fective. In the small temperature difference limit, the 
Prandtl number of the modeled fiuid can vary in a rela- 
tively large range, comparing with the work of Bouchut 



and Perthame 



the continuous BGK model whose 



Prandtl number may vary only from to 1. Applica- 
tion of this lattice BGK model to the more complicated 
Rayleigh-Benard fiow is under current research work. 
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Fig. 1. Comparison of the numerical measured Prandtl 
numbers with the analytical prediction with formula (7). 
"2D16VP1" and "2D16VP2" are code names for different 
models. "2D16V" stands for two dimensions and sixteen 
velocity model, while "PI" represent a plain cr-incoporated 
model and "P2" a model with a further modification to min- 
imize the deviation. 
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Fig. 2. Temperature pofiles obtained from lattice BGK 
simulation of the Couette flow. The lattice BGK model used 
here has 17 velocities, including the stationary particle dis- 
tribution. 
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